﻿using System;

namespace TribesPSO
{
    /// <summary>
    /// A random number generator implementation capable of generating random numbers based
    /// on a Gaussian distribution
    /// </summary>
    public class GaussianRandom : Random, IGaussianRandom
    {
        private const double c_0 = 2.515517;
        private const double c_1 = 0.802853;
        private const double c_2 = 0.010328;
        private const double d_1 = 1.432788;
        private const double d_2 = 0.189269;
        private const double d_3 = 0.001308;

        /// <summary>
        /// Initializes a new GaussianRandom pseudo random number generator object seeded with the system time
        /// </summary>
        public GaussianRandom() : base() { }

        /// <summary>
        /// Initializes a new GaussianRandom pseudo random number generator object seeded with the specified value
        /// </summary>
        /// <param name="seed">A number used to calculate a starting value for the pseudo-random number sequence</param>
        public GaussianRandom(int seed) : base(seed) { }

        /// <summary>
        /// Choose the next random number.  The distribution of randomness is "normal"
        /// or "Gaussian" with a mean value 0 and a standard deviation of 1
        /// </summary>
        /// <returns>A random value based on the gaussian distribution</returns>
        public double NextGaussian()
        {
            return NextGaussian(0, 1);
        }

        /// <summary> 
        /// Choose the next random number. The distribution of randomness is "normal" 
        /// or "Gaussian" with a mean value dMu and a standard deviation of dSigma 
        /// </summary> 
        /// <param name="mu">The center point (mean) of the normal distribution</param> 
        /// <param name="sigma">The standard deviation of the distribution</param> 
        /// <returns>Random value</returns> 
        public double NextGaussian(double mu, double sigma)
        {
           /* As of 12/2010 this function is the slow path in the PSO algorithm (excluding the objective function
            * for certain problems).  You could probably improve performance by replacing this gaussian function
            * with the zigurrat algorithm.  Remember to profile changes so you don't accidentally slow things down!
            */
            return mu + CumulativeGaussian(base.NextDouble()) * sigma;
        }

        /*There's a nasty bug here.  Although it's unlikely, random.NextDouble() can return 0.  In that case, p will be 0
         * and when we get to the line:
         * double t = Math.Sqrt(Math.Log(1.0 / (p * p)));
         * we're going to get a suprise!
         * The upside is that when I say it's unlikely, I mean SUPER unlikely.
         */
        private static double CumulativeGaussian(double p)
        {
            // p is a rectangular probability between 0 and 1 
            // convert that into a gaussian. 
            // Apply the inverse cumulative gaussian distribution function 
            // This is an approximation by Abramowitz and Stegun; Press, et al. 
            // See http://www.pitt.edu/~wpilib/statfaq/gaussfaq.html 
            // Because of the symmetry of the normal 
            // distribution, we need only consider 0 < p < 0.5. If you have p > 0.5, 
            // then apply the algorithm below to q = 1-p, and then negate the value 
            // for X obtained. 
            bool fNegate = false;

            if (p > 0.5)
            {
                p = 1.0 - p;
                fNegate = true;
            }

            double t = Math.Sqrt(Math.Log(1.0 / (p * p)));
            double tt = t * t;
            double ttt = tt * t;
            double X = t - ((c_0 + c_1 * t + c_2 * tt) / (1 + d_1 * t + d_2 * tt + d_3 * ttt));
            return fNegate ? -X : X;
        }

    }

}